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Abstract 

This article proposes spectral numerical methods to solve the time evolution of a convection problem with 
viscosity depending exponentially on temperature. The set-up is a 2D domain with periodic boundary conditions 
along the horizontal coordinate. At a fixed aspect ratio, the analysis is assisted by bifurcation techniques such 
as branch continuation, which has proven to be a useful, systematic method for gaining insight into the possible 
stationary solutions satisfied by the basic equations. Stable stationary solutions become unstable through a 
Hopf bifurcation, after which the time-dependent regime is solved by the spectral techniques proposed for that 
purpose in this article. The morphology of the plume is described and compared with others obtained in the 
literature. 
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1 Introduction 

Thermal convection of the upper mantle plays an important role in Earth dynamics [H [2j 13 , 3], since observations 
indicate that on geologic time scale the mantle deforms like a fluid. If the mantle is treated as such, properties 
such as thermal conductivity, density or viscosity v must be modelled. For instance, in the upper mantle viscosity 
has typically been considered to depend exponentially on temperature [SJ |5J [7] or to follow the Arrhenius law [S] . 
Interest in the numerical study of these problems is clear for its geophysical applications. On the one hand from 
the mathematical point of view it has recently been proven that this is a well posed problem j9[ [10] for dependences 
which are smooth bounded positive analytical functions, so it stands a good chance of solution on a computer 
using a stable algorithm. The variability in viscosity introduces strong couplings between the momentum and heat 
equations, as well as introducing important nonlinearities into the whole problem. On the other hand, it is of major 
interest in mantle convection to consider the fact that the Prandtl number, which is the quotient of viscosity and 
thermal diffusivity, is virtually infinite. This limit transforms the set of equations describing the time-dependent 
problem into a differential algebraic problem (DAE), which is very stiff. 

In this context, this article discusses the performance of several time evolution spectral schemes for a convec- 
tion problem with viscosity dependent on temperature. In particular, we extend the results discussed in [7J by 
considering an exponential law in a 2D domain with periodic boundary conditions along the horizontal coordinate. 
We characterize time-dependent solutions demonstrating the efficiency of the time-dependent scheme to describe 
the solutions beyond the stationary regime. 

Spectral schemes are widely used to numerically integrate solutions of convection problems with constant 
viscosity QTJ [T21 EH [H]. They are preferred to other discretization methods because they are highly accurate 
and rapidly convergent 115 . In thermoconvection problems, pressure is usually avoided by using, for instance, 
potentials of velocity [T31 [16l [17] . This methodology eliminates pressure from the equations and replaces the 
velocity with scalar fields that raise the order of the differential equations, thus requiring additional boundary 
conditions. Alternatively, solving thermoconvection in the primitive variables formulation introduces the problem 
of handling first order derivatives for pressure in a closed container and determining its boundary conditions. In 
the velocity-pressure formulation [THl Q33 HD] > first order derivatives are avoided by deriving a Poisson equation for 
pressure |18) . although the problem of finding its boundary conditions still remains. An efficient spectral method 
for the time evolution of convection problems with constant viscosity in primitive variables is presented in [12]. A 
different and also successful spectral method in primitive variable formulation is discussed in [211 G2] [23] , but only 
for stationary problems. 

Spectral methods in the mantle convection community are not so popular [53], as they are reported to have 
limitations when handling lateral variations in viscosity. Alternatively, preferred schemes exist in which the basis 
functions are local; for example, finite difference, finite element and finite volume methods. For instance, the works 
by [25] [26] have treated this problem in a finite element discretization in primitive variables, while in [27j [28] 



finite differences or finite elements are used in the stream-function vorticity approach. Spectral methods have been 
successfully applied to model mantle convection with moderate viscosity variations in, for instance, [351 [30]. These 
works do not use the primitive variables formulation, and deal with the variations in viscosity by decomposing it 
into a mean (horizontally averaged) part and a fluctuating (laterally varying) part. Our approach addresses the 
variable viscosity problem by proposing a spectral approximation in primitive variables without any decomposition 
on the viscosity. The main novelty in this paper is the extension of the spectral methodology discussed in [7] , valid 
only for stationary problems, for solving the time-dependent problem, and the extension of the results to describe 
time-dependent solutions. 

The use of spectral methods in mantle convection as an alternative to local basis of functions is specially 
interesting for characterizing solutions in the presence of symmetries. This is of particular interest in our problem, 
where the equations with periodic boundary conditions are invariant under horizontal translations, represented by 
the SO(2) group. Assemat and co-authors [3T] have used high order finite element methods to solve a constant 
viscosity convection problem in a circular container. These authors have noted the influence of the computational 
grid on the breaking of the 0(2) symmetry, by producing pinning effects on the solution. They have also noted that 
heteroclinic cycles reported in similar set-ups, both experimentally [32] and with a fully spectral approach [33], are 
solutions absent in their computational approach. Heteroclinic cycles are known to be persistent in systems with 
0(2) symmetry as discussed in [Ml [351 • Thus, a motivation for the use of a spectral approach is that it would 
potentially allow the description of solutions which may be overlooked with approaches based on finite elements. 




Figure 1: Problem set-up. 



As regards temporal discretization, backward differentiation formulas (BDF's) are widely used in convection 
problems. This is the case of the work discussed in [12], which following ideas proposed in [37] uses a fixed time 
step second-order-accurate which combines Adams-Bashforth and BDF schemes. A recent article by Garcia and 
co-authors t I3] compares the performance of several semi-implicit and implicit time integrations methods based on 
BDF and extrapolation formulas. The physical set-ups discussed in these papers are for convection problems with 
constant viscosity and finite Prandtl number. In contrast, this article focuses on convection problems with viscosity 
strongly dependent on temperature and infinite Prandtl number that lead to a differential algebraic problem. We 
will see that the semi- implicit methods discussed in [371Q2] do not work in this context. BDFs and implicit methods 
are known to be an appropriate choice [351 EHj for efficiently tackling very stiff problems. According to [351 113) . 
for the time discretization scheme we propose several high order backward differentiation formulas which are ready 
for an automatic stepsize adjustment. Furthermore, we solve the fully implicit problem and also propose a semi- 
implicit approach. The output and performance of this option are compared with those of the implicit scheme. It 
is found that the semi-implicit approach presents some advantages in terms of computational performance. 

The article is organized as follows: In Section[2j we formulate the problem, providing a description of the physical 
set-up, the basic equations and the boundary conditions. Section [3j describes a spectral scheme for stationary 
solutions which will be useful for benchmarking the time dependent numerical schemes. First the conductive 
solution and its stability is determined. Other stationary solutions appear above the instability threshold, which 
are computed by means of a Newton-Raphson method using a collocation method. The stability of the stationary 
solutions is predicted by means of a linear stability analysis and is also solved with the spectral technique. Section 
[4] discusses several time-dependent schemes, which include implicit and semi-implicit schemes. Section [5] reports 
the results at a fixed aspect ratio in a range of Rayleigh numbers. Stationary and time-dependent solutions are 
found and different morphologies of the thermal plumes are described. Some computational advantages of some 
schemes versus others are discussed. Finally, Section [6] presents the conclusions. 
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2 Formulation of the problem 



The physical set-up, shown in Fig[T] consists of a two dimensional fluid layer of depth d (z coordinate) placed 
between two parallel plates of length L. The bottom plate is at temperature To and the upper plate is at T±, where 
T\ = To — AT and AT is the vertical temperature difference, which is positive, i.e, T± < T). 

In the equations governing the system, u = (u Xl u z ) is the velocity field, T is the temperature, P is the pressure, 
x and z are the spatial coordinates and t is the time. Equations are simplified by taking into account the Boussinesq 
approximation, where the density p is considered constant everywhere except in the external forcing term, where a 
dependence on temperature is assumed as follows p = po{\ — a(T — T±)). Here po is the mean density at temperature 
Ti and a the thermal expansion coefficient. 

We express the equations with magnitudes in dimensionless form after rescaling as follows: (x' , z') — (x, z)/d, 
f = Kt/d 2 , u' = da/K, P' = d 2 P/(p Kis ) , 9' = (T — Tl)/(AT). Here k is the thermal diffusivity and v 
is the maximum viscosity of the fluid, which is the viscosity at temperature T\. After rescaling the domain, 
fii = [0, L) x [0,d] is transformed into = [0, T) x [0, 1], where T — L/d is the aspect ratio. The system evolves 
according to the momentum and the mass balance equations, as well as to the energy conservation principle. The 
non-dimensional equations are (after dropping the primes in the fields): 

V-u = 0, (1) 

/a., i V7..\ _ r>o^> nn i j:.. I \ ) /w.. i lx-r..\l 1 



■(ftu + u.Vu) = i?0e 3 - VP + div -^(Vu+(Vu) J ) , (2) 
Pr \ v Q J 

d t 9 + u • V0 = A0. (3) 

Here e*3 represents the unitary vector in the vertical direction, R = d 3 agAT '/ (vqk) is the Rayleigh number, g is the 
gravity acceleration and Pr = vq/k is the Prandtl number. Typically for rocks, Pr is very large, as they present 
low thermal diffusivity (approximately 10~ 6 m 2 /s) and very large viscosity (of the order 10 20 Ns/m 2 ) [5]. Thus, 
for the problem under consideration, Pr may be considered as infinite and the term on the left-hand side in ([2| 
can be made equal to zero. The viscosity v(0) is a smooth positive bounded function of 9. For describing mantle 
convection problems, a usual choice for v(9) is an exponential law El [7]; according to the results discussed in 
[7j, this is used in dimensionless form: 

^ = cxp(-pR9) (4) 

Dependence on Eq. Q reduces to that of constant viscosity if p, — 0, while for temperature-dependent viscosity 
we consider p = 0.0862. The viscosity contrast in equation Q for Rayleigh numbers up to R = 120 -as employed 
in this article- is 3.1 • 10 4 . 

For boundary conditions, we consider that the bottom plate is rigid and that the upper surface is non deformable 
and free slip. The dimensionless boundary conditions are expressed as, 

9 = 1, u = 0, on z — and 9 = d z u x — u z — 0, on z = 1. (5) 

Lateral boundary conditions are periodic. Jointly with equations ([T])-([3|, these conditions are invariant under 
translations along the ^-coordinate, which introduces the symmetry SO(2) into the problem. In convection problems 
with constant viscosity, the reflexion symmetry x — > —x is also present insofar as the fields are conveniently 
transformed as follows (9,u x ,u z ,p) — > (9, —u x ,u z ,p). In this case, the 0(2) group expresses the full problem 
symmetry. The new terms introduced by the temperature dependent viscosity, in the current set-up Eq. ^ 
maintain the reflexion symmetry, and the symmetry group is 0(2). 



3 Numerical schemes for stationary solutions 

The numerical codes proposed in this article for the time-dependent problem requiere a priori known solutions 
for benchmark. In the set-up under consideration, both stationary and time-dependent solutions exist. Some of 
the stationary solutions have simple analytical expressions which are known a priori, but there are others which 
are not so simple and must be found numerically. Stationary solutions are not stable in the full parameter space, 
and after becoming unstable either other stationary solutions may become stable or a time-dependent regime is 
observed. In order to verify our methods, the self-consistency of results provided by a different kind of analysis is 
required. In this context, this section describes stationary solutions to the system ([l])-([3]) and their stability. This 
description, together with the fact that the problem is well-posed |10j , provide a required a priori knowledge that 
will assist in the examination of the validity of different time-dependent numerical schemes. 
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3.1 The conductive solution 



The simplest stationary solution to the problem described by equations ([lJ)-((3]) with boundary conditions ([5]) is 
the conductive solution which satisfies u c = and 9 C = — z + 1. This solution is stable only for a range of vertical 
temperature gradients which are represented by small enough Rayleigh numbers. Beyond the critical threshold R c , 
a convective motion settles in and new structures are observed which may be either time dependent or stationary. 
The stationary equations in the latter case, obtained by canceling the time derivatives in the system Q-Q, are 
satisfied by the bifurcating solutions. As in the conductive solution the new solutions depend on the external 
physical parameters, and new critical thresholds exist at which they lose their stability, thereby giving rise to new 
bifurcated structures. 




r 



Figure 2: Critical instability curves R(m,T) for a fluid layer with temperature dependent viscosity /j, = 0.0862 
according to the exponential law given in Q . 

In this section, we first describe the instability thresholds for the conductive solution. For this purpose small 
perturbations are added to it: 



u(x, z, t) 
9(x, z, t) 
P{x, z, t) 



+ u(z)e xt+ikx , 

-z + l + 9(z)e xt+tkx , 

-Rz 2 /2 + Rz + C + P{z)e xt+lkx . 



(6) 
(7) 
(8) 



The sign in the real part of the eigenvalue A determines the stability of the solution: if it is negative the perturbation 
decays and the stationary solution is stable, while if it is positive the perturbation grows in time and the conductive 
solution is unstable. If these expressions are introduced into the system Q-([3]), and both the nonlinear terms in 
the perturbations and their tildas are dropped, the system becomes: 








xe 



iku x + d z u 
ikP- ' 6 > V M 

d z P 

(C - 



"0 

k 2 )9 + u z 



(iku z - 
-d z u z 



d z u x ) + ^(d 2 zz 



^0 



(d\ z - k 2 )u z 



k 2 )u x 



R6, 



The boundary conditions for the perturbation fields are: 

9 = 0, u = 0, on z = and ( 



= d z u x = u z = 0, on 2 = 1. 



(9) 



Stability analysis of the conductive solution for viscosity dependent on temperature is numerically addressed in 
[7j |40] , although for the results reported below we follow the scheme presented in [21] . Appropriate expansions of 
the unknown fields (u, 9, P) transform the eigenvalue problem into its discrete form: 

Aw = XBw, 



where the expansion coefficients are stored in the vector w. This generalized eigenvalue problem supplies the 
dispersion relation R = R(k) at the bifurcation point (Re(A) = 0). The wavenumber k in this relation cannot be 
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arbitrary but must meet the periodic boundary conditions, i. e. 



kT = 2irm, with m = 1, 2, 3 . . . (10) 
Here, to is the number of wavelengths of the unstable structure growing in the finite domain. Restrictions to k 



given by condition (10 1 are replaced in the dispersion relation R = R(k), thereby providing a critical curve for each 
integer to as a function of the aspect ratio r, R — R(m,T). Figure [2] displays critical Rayleigh numbers R c on the 
vertical axis as a function of the aspect ratio T on the horizontal axis. The conductive solution is stable below the 
critical curves, which means that in a box with a given aspect ratio, r, if R < R c then initial conditions near to 
the conductive solution evolve in time approaching it. Alternatively, if at that aspect ratio R > R c , then initial 
conditions which are near to the conductive solution evolve in time away from it, towards a different solution. This 
new solution may be stationary or time-dependent. Figure [2] confirms that for increasing aspect ratios the most 
unstable spatial eigenfunction increases its wavenumbcr m. 

The critical Rayleigh number in the problem under study is very low (R c ~ 74) if compared with that obtained at 
constant viscosity or with that reported in problems in which viscosity depends on temperature. The reason for this 
is that the viscosity isq used to define the Rayleigh number is the maximum viscosity at the upper boundary. This 
viscosity is between two to four order of magnitudes larger than the lower viscosity in the fluid, which furthermore 
is the viscosity typically used to define the R number in other works. For this reason, our critical R number may 
be between two to four order of magnitudes smaller than that in previously reported results. 



3.2 Numerical stationary solutions 

There exist stationary solutions to the system ([l])-([3]) that bifurcate from the conductive solution above the insta- 
bility thresholds displayed in Figure 2 and may be computed numerically. They are stationary because they satisfy 
the stationary version of equations ( 1 )-([3]), which are obtained by canceling the partial derivatives with respect to 
time. These solutions may be numerically obtained by using a variant of the iterative Newton-Raphson method, 
similar to that in [7]. This method starts with an approximate solution at step s = 0, to which is added a small 
correction in tilda: 

(ir 5 + u,9 s + 9,P S + P). (11) 

These expressions are introduced into the system |l])-([3|, and after canceling the nonlinear terms in tilda, the 
following equations are obtained: 



=V • u + V • u s , 

= - d x P - d x P s + — [L n (d s ,K,ul) + L 12 (9 s )u x + L 13 (8 s )u z + L 14 (9 S , u%, u s z )6], 

= -8 z P- d z P s + - [L 21 (9 S , u%, u%) + L 22 (8 s )u x + L 23 (9 s )u z + (L 24 (9 S ,u%,u s z ) + R)0], 
=u • V(9 S + u s • V9 + u s ■ V6» s -A9- A9 S . 



(12) 
(13) 

(14) 
(15) 



Here, Lij (i — 1, 2, j — 1, 2, 3, 4) are linear operators with non constant coefficients which are defined as follows: 



L n (( 


),u x ,u z ) 


=2d e v{9)d x 9d x u x 


+ v(9)Au x + d e v(9)d z 9{d x u z + 


d z u x ), 




(16) 




Ly 2 {8) 


=2d e v{9)d x 9d x + 


v(9)A + d e v(0)d z 6d x , 






(17) 




Lu(9) 


=d e v{6)d z 6d x , 








(18) 




l,u x ,u z ) 


=2d s v(9)d x u x d x - 


f 2d 2 09 u{9)d x 8d x u x + d e v(8)Au x 


+ {d x u z + d z u x )(dev{8)d z - 


h d 2 ee v{9)d z 9), 


(19) 


L 21 (t 


)>u x ,u z ) 


=2d e v{9)d z 9d z u z 


+ v{6)Au z + d e u{9)d x 9{d z u x + 


d x u z ), 




(20) 




L 22 {9) 


=d e v{6)d x 6d z , 








(21) 


£23 (< 


l,u x ,u z ) 


=2d e v{9)d z 9d z + 


v{0)A + d e v{6)d x 0d z , 






(22) 


L 24 (t 




=2d e v(9)d z u z d z - 


V 2d 9e v{6)d z 6d z u z + d e u{9)Au z 


+ (d z u x + d x u z )(dgv(9)d x 4 


d e ev(6)d x 6). 


(23) 



The unknown fields u, P, 9 are found by solving the linear system with the boundary conditions ([9|) and the new 
approximate solution s + 1 is set to 



u. 9 



s+1 



6, P 



s+1 



ps 



The whole procedure is repeated for s + 1 until a convergence criterion is fulfilled. In particular, we consider that 
the I 2 norm of the computed perturbation should be less than 10 -9 . 
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At each step, the resulting linear system is solved by expanding any unknown perturbation field Y, in Chebyshev 
polynomials in the vertical direction and Fourier modes along the horizontal axis: 



Y(x,z) 



\L/Z] M-l 

E E 

1=1 m=0 



L 



M-l 



l=\L/2\+l m=0 



aim T m{z)e 



i(l-l-L)x 



(24) 



In this notation, [•] represents the nearest integer towards infinity. Here L and M are the number of nodes in 
the horizontal and vertical directions, respectively. Chebyshev polynomials are defined in the interval [—1, 1] and 
Fourier modes in the interval [0, 2ir). Therefore, for computational convenience, the domain f2 = [0,r) x [0,1] is 
transformed into [0, 27r) x [— 1, 1]. This change in coordinates introduces scaling factors into equations and boundary 
conditions which are not explicitly given here. There are 4xLxM unknown coefficients which arc determined 



by a collocation method in which equations ( 12 1-( 14 ) and boundary conditions are posed at the collocation points 



Uniform grid: Xj = (j — 1) 
Gauss-Lobatto: z. 



cos 



2tt 

T' 

i - 1 
M-l 



-If 



j = 1, . . . ,L; 

i = l,...,M; 



After replacing expression (24 1 in equations ( 12 )- ( 14 1 , the partial derivatives are evaluated on the basis of functions. 



Derivatives of the Chebyshev polynomials at the collocation points are evaluated a priori according to the ideas 
reported in |41j . The expansion (24) is an interpolator outside the collocation points, which when restricted to 
these points may be rewritten as: 



L M-l 

E E 

1 = 1 m=0 



(25) 



for I > L/2 at the collocation points. However, 



due to the aliasing effect of the functions e l ^ L ~ 1 ' >Xj and e *('~ 1_ - £ ') x J 
this expression is not valid for computing the spatial derivatives of the fields, for which purpose expansion ( 24 ) 
should be used. Although in practice Eq. (24) is correct and provides good results, we do not employ it in this 



work because it involves complex functions and complex unknowns af m , eventually leading to the inversion of 
complex matrices which computationally are more costly than real matrices. On the other hand since the unknown 
functions Y are real, they admit expansions with real functions and real unknowns. In order to obtain these 
functions, we take into account Euler's formula e %lx = cos(Zx) + isin(lx), which is replaced in (24). We also note 
that the coefficients are given by conjugated pairs in such a way that, for instance, a\ m = a^* t . Strictly speaking, 



expansion (24) is a real function only if every coefficient in the first summatory for I > 2 has a conjugate pair in 



the second summatory. This implies that L must be an odd number; thus in what follows we restrict ourselves to 
odd L values. With these considerations the following equations are obtained: 



Y{x,z) 



[L/2] M-l 

E E bf m T m {z)cos{{l-l)x) 



1=1 m=0 



[L/2] M-l 

E C Ym T m(z) Sin((l - l)x). 
1=2 m=0 



E 



(26) 



ij m and bj m = 23?(aL) and cl 



-29f(o 



Some relations among the real and complex coefficients are: b\ m 
for 1 = 2,..., [L/2]. 

The rules followed to obtain as many equations as unknowns are described next. Equations ( 12 )— ( 15 ) are 
evaluated at nodes i = 2, ...,M — 1, j = 1,...,L. This provides 4 x (M — 2) x L equations; the boundary 
conditions ([5| are evaluated at i = 1, M, j = 1, . . . , L. This supplies additional 4 x L x M — 2L equations. In order 
to obtain the remaining 2L equations, we complete the system with extra boundary conditions which eliminate 
spurious modes for pressure |21l I42j projecting the equation of motion in the upper and lower plate of the domain 
i.e. the equation (14) is evaluated at nodes i = 1, M, j = 1, . . . , L. This choice has been reported to be successful 



for many convection problems [21, 151 144j . However, in the present set-up, results are improved if, for the equation 



(14) imposed at the upper boundary, the continuity equation is assumed and d Z7 u z is replaced with —d xz u x . With 
these rules we obtain a linear system of the form AX = b in which X contains the unknowns. However, the matrix 
A is singular due to the fact that pressure with the imposed conditions is defined up to an additive constant. 
According to[3Tl 221 |U), we fix the constant by removing equation (14) at node j = 1, i = 2 and adding at this 
point the equation bf = 0. This is computationally cheaper than the pseudo-inverse method. However, in the 



problem under study, dropping the equation ( 14 ) at one point introduces weakly oscillating structures on this side 
of the pressure field. To overcome this drawback, in the final step of the iterative procedure, once the tolerance 
is attained, we proceed alternatively by computing a pseudoinverse of the matrix A by using the singular value 
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decomposition (SVD). Let A = UT.V* be the singular value decomposition of A; VS + U* is the pseudoinverse of 
A, where S + is the pseudoinverse of a diagonal matrix, i.e. it takes the reciprocal of each non-zero element on the 
diagonal, and transposes the resulting matrix. 

The study of the stability of the numerical stationary solutions under consideration is addressed by means of a 
linear stability analysis. Now perturbations are added to a general stationary solution, labeled with superindex b: 



u(x,z,t) = u b (x, z) + u(x, z)e xt , 
9(x,z,t) = 9 b {x,z) + 9(x,z)e xt , 
P(x,z,t) = P b (x,z) + P(x,z)e xt . 



The linearized equations are: 

=V ■ u 

= - d x P 

= -d z P- 

o =u • ve b 



~[L 12 (9 b )u x 

~[L 22 (9»)u x 
u b ■V9 + u b ■ 



+ L 23 {9 b , 
V9 b - A9 



L u {9 b ,u b x ,u b z ) 



u, 



+ (£ 24 ( 
X9, 



R)9] 



(27) 
(28) 
(29) 



(30) 
(31) 

(32) 
(33) 



where the operators Lij are the same as those defined in Eqs. (16 1- (23). The stability of the stationary solutions is 



approached with the collocation method used for the Newton-Raphson iterative method. Expansions of the fields 
(26) are replaced in equations (30)-(33), and they and the boundary conditions Q are evaluated at the collocation 
nodes following the same rules as before. As a result, the discrete form of the generalized eigenvalue problem is 
obtained: 



Aw = XBw, 



(34) 



where w is a vector containing unknowns. 

In order to solve the generalized eigenvalue problem, we use a generalized Arnoldi method, which is described 
in [IS]. The numerical approach uses the idea of preconditioning the eigenvalue problem with a modified Caley 



transformation, which transforms the problem (34) and which admits infinite eigenvalues into another one with all 



its eigenvalues finite. Afterwards, the Arnoldi method is applied. 



4 Numerical schemes for time dependent solutions 

The governing equations ([T])-([3]), together with boundary conditions ([5]), define a time-dependent problem for which 
we discuss temporal schemes based on a primitive variables formulation. The spatial discretization is analogous 
to that proposed in the previous section, thus the focus in this section is to discuss the time discretization of the 
problem. 

To integrate in time, we use a third order multistep scheme. In particular, we use a backward differentiation 
formula (BDF), since as discussed in 38, 3S] these are highly appropriate for very stiff problems such as ours. The 
BDF evaluates the time derivative in Eq. ([3| by differentiating the formula that extrapolates the field 9 n+1 with 
a third order Lagrange polynomial that uses fields at times 9 n ,9 n ~ 1 ,9 n ~ 2 , i.e.: 

9{t) := e n+1 (t)9 n+1 + £ n (t)9 n + in-iit)^- 1 + e n - 2 (t)9 n - 2 

Ik being a Lagrange polinomial of order 3: 



(i)= [] n-2<k<n + l 



n — 2 < to < n 
m 7^ k 



Then, 

d t 9 n+1 = i' n+1 {t n+x )6 n+1 +e' n (tn+i)9 n +l' n _ l {t n+1 )6 n - 1 +£' n _ 2 (t n+1 )9 n - 2 (35) 
For the fixed time step case, the time differentiation simplifies to the equation: 

x = iig»+ 1 -i8g» + 9e»- 1 -2g»- 2 
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Figure 3: Bifurcation diagram for a fluid with viscosity dependent exponentially on temperature (fi = 0.0862 in 
law Q) at r = 3.4. Stable branches are solid while unstable branches are dashed. Solutions are displayed at the 
R numbers highlighted with vertical dashed lines. 



In stiff problems, a variable time step scheme with an adaptative step control that adjusts itself conveniently to 
the different regimes is advisable. The expressions for the time derivative (35) for the variable time step case are: 

At 2 n + 4At n+1 At n + At„At„_i + 3A4 +1 + 2Ai n _ 1 At n+1 



^1-2(^+1) 



(Ai„ + At n+1 + At n _!)(A£„ + At n+1 )At n+1 
Atj + A£„At ra _i + At n -iAt n+ i + Atl +1 + 2At n At n+1 

(At n + At n ^)At n At n+1 
At n+1 (At n + At n+1 + At n _i) 
(At n + At n+ i)At n At n -i 

. At n+1 (At n + At n+1 ) 

" [(Atn + Ai„_i) 2 + At n At n+1 + At„_!At n+1 )Ai„_i ' 



where A£„ = £„ — t n —\. 

The variable time step scheme controls the step size according to the general ideas proposed by j46l |38] , and 
adapted to our particular case with parameters taken from [47) . The result of an integration at time n + I is 
accepted, depending on the estimated error E for the fields. The error estimation E is based on the difference 
between the solution obtained with a third and a second order scheme. Then essentially the new time step is 
evaluated as follows: 

/ E \ -1/(9+1) 

h n ew = s — i hold- 

\ tolerance J 

In practice, this expression is tuned and a maximum increase of the step size is allowed. Here s is a safety factor, 
and q is the order of the numerical scheme. Acceptance of the result of an integration means that E is below a 
certain tolerance, which is explained in the subsection |4.2| If the result is accepted, a new time step is proposed 
according to the law: 

h I «holi (tou^Y . E > tolerance ■ (f ) 1/p 

{ 5h otd , E< tolerance - (§) , 

where h is the size of the time step and s = 0.9, p = —0.33 (q = 2). In case of rejection, the step is decreased as 
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follows: 



with q = 3. 



\ tolerance / 



-1/(9+1) 



old • 



4.1 The fully implicit method 

BDFs are a particular case of multistep formulas which are implicit, thus the BDF scheme implies solving at each 
time step the problem (see |38j): 

d t y n+1 = f(y' l+1 ), 



which in the particular problem under consideration becomes: 



= V ■ u 



n.+ l 



= R6 n+L e 3 - VP 



71+ 1 



div 



qn+l 



(Vu" +i + (Vu n+ f) 



'A) 



Q td n+1 = _ u n+l . v6 p+l + A6 ,«+l ; 



(37) 
(38) 

(39) 



where dtO n+ is replaced by the expression (35). The solution to the system ( 37 )-( 39 ) is our benchmark for 
transitory and time- dependent regimes. 

The nonlinear terms on the right-hand side of these equations are approached at each time step, n + 1, by a 
Newton-Raphson method similar to the one described before to find numerical stationary solutions. We assume 
that the solution at time n + 1 is a small perturbation Y of an approximate solution. Linear equations for Y 



are derived by introducing the analogue of expression (111 into the nonlinear terms of equations ( 37 1-( 39 ) and 



cancelling all the nonlinear terms in tilda. The resulting linearized terms are the same as those appearing in Eqs. 
(12|-(15|. For the first step s = 0, we ta ke a s an approximate solution the solution at time n. The unknown 
perturbation fields Y are expanded by Eq. (26), which is replaced in the equations and boundary conditions at the 
collocation points according to the rules described in the previous section. In particular, the improved boundary 
conditions for pressure, previously described, are used. A linear system is obtained: 



Ay 



(40) 



which is iteratively solved at each time step n + 1 until the perturbation Y is below a tolerance. Here, A is a 
matrix of order 4 x L x M and y is the vector containing the unknown coefficients of the fields Y. As previously 



performed the matrix A is converted into a full rank matrix by removing the projection of equation (38) at node 
j = 1, i = 2, and fixing the pressure by adding at this point the equation bf Q = 0. This provides accurate results 
and is computationally cheaper than the pseudoinverse method. 



4.2 The semi-implicit method 

Implicit methods are a robust and numerically stable choice for stiff problems. However, they may be rather 



demanding computationally, since at each time step they require several matrix inversions to solve the system ( 40 ) 
at successive iterations. This subsection proposes an alternative semi-implicit scheme that is computationally less 
demanding than the fully implicit scheme. 



Similarly to the fully implicit method, the semi-implicit scheme approaches the nonlinear terms in Eqs. ( 37 )-( 39 ) 
by assuming that the solution at time n + 1 is a small perturbation Z of the solution at time n; thus, z n+1 = z n + Z. 
Once linear equations for Z are derived, the equations are rewritten by replacing Z = z n+1 — z™. Additionally, 



the linear system is completed by using expression (35 1 for the time derivative of the temperature. The solution is 



obtained at each step by solving the resulting linear equation for variables in time n + 1. As before, the unknown 



fields are expanded by Eq. (26), which is replaced in the equations and boundary conditions at the collocation 



points according to the rules described in the previous section. At each time step the linear system: 

Ay = b (41) 

is solved, where A is a matrix of order 4xLxM and y is the vector containing the unknown coefficients of the 
expansions of the z n+1 fields. The matrix is transformed into a full rank matrix with the same procedure used for 
the fully implicit case. 

The fully implicit and the semi-implicit methods are implemented with a variable time step scheme that requires 
an error estimation based on the difference between the solution obtained with a third and a second order scheme. 
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M=36 


M=38 


M=40 


M=42 


L=29 


0.0030 - 0.0002i 


0.0046 


0.0046 


0.0047 




-8.2205 


-8.4419 


-8.4419 


-8.4419 


L=31 


0.0017 


-0.0019 


0.0017 


0.0017 




-8.4418 


-8.4417 


-8.4418 


-8.4418 


L=33 


0.0011 


0.0011 


-5.79o-4 -l.llc-4i 


7.95c-4 + 6.30o-5i 




-8.4418 


-8.4418 


-8.4417 


-8.4418 



Table 1: Computation of the two eigenvalues with largest real part for the stationary solution obtained at T = 3.4, 
R = 78 at different expansions L x M. 





M=30 


M=40 


M=50 


M=60 


M=70 


M=80 


L- 


=31 


0.9433 + 0.7023i 


0.2227 


0.5829 + 0.4345i 


0.3926 


0.4217 


0.2265 






0.9433 - 0.7023i 


-2.9525 


0.5829 - 0.4345i 


-2.3225 + 0.7449i 


-2.1744 + 0.7450i 


-2.6555 + 0.7536i 


L 


=37 


0.3564 


0.1837 


0.1369 


0.1318 


0.1834 


0.0946 






-2.1572 + 0.0829i 


-2.1273 + 0.0843i 


-2.6814 


-2.0368 + 0.0843i 


-2.1427 + 0.0841i 


-2.1495+ 0.0964i 


L 


=43 


-0.2039 + 0.0913i 


0.1767 


0.0773 +0.037H 


-0.0289 + 0.0886i 


0.0585 


0.0464 






-0.2039 - 0.0913i 


-2.7849 + 0.0546i 


0.0773 -0.0371i 


-0.0289 - 0.0886i 


-2. 1576+9. 7059e-3i 


-2.1308 


L 


=49 


-0.1436 


0.0836 


0.0705 


0.0379 


9.7501e-3 


4.7754e-3 






-2.4633 +0.8885i 


-2.8147 + 0.0127i 


-2.3698 


-2.3625 


-2.1391+ 6.3240e-4i 


-2.1482 


L- 


=55 


-0.1382 


0.0571 


0.0315 


5.8407e-3 


3.7262e-3 


7.8574e-4 






-2.5619 + 0.8760i 


-2.1569 


-2.1893 + 1.0669e-3i 


-2.1618 + 6.7549e-4i 


-2.1597 


-2.1598 


L- 


=61 


0.0913 


0.0261 


4.5781e-3 


1.5021e-3 


6.4589e-4 


9.8951e-5 






-2.6323 + 0.1950i 


-2.1623 


-2.1606 


-2.1603 


-2.1599 


-2.1599 



Table 2: Computation of the two eigenvalues with largest real part for the stationary solution obtained at T = 3.4, 
R = 110 at different expansions L x M. 



The second order scheme approaches the time derivative of the temperature as follows: 

d 9 „+i = (Ag + 2At n At n+1 )8 n + 1 - (At n + At ra+1 ) 2 0» + Atl +1 0^ 
* ~ At n At n+1 (At n+1 + At n ) ' 1 ' 

where At n = t n — f n _i. In practice, this changes the linear system to be solved at each step, as follows: 

Ay = b 

The computation of the second order solution thus leads to an additional matrix inversion at each time step, and 
as a consequence the full calculation is slowed down considerably. To avoid this additional inversion, we estimate 
the error by measuring instead how well the third order solution y satisfies the second order system, i.e.: 

_ \\b-Ay\\ 

E " (43) 

where || • |j represents the l 2 norm. Acceptance of the result of an integration means that E is below a tolerance 
that we fix at 5 • 10~ 6 . Once the error is estimated, the step size is determined as explained at the beginning of 
Section 4. 



4.3 Other semi-implicit schemes 

For completeness, we describe here alternative semi-implicit schemes that have been successfully used in fluid 
mechanics and convection problems with constant viscosity. Despite its success in many fluid mechanics set-ups, 
this scheme is insufficient for our problem. 

The adaptation of the numerical scheme described in [371 112j to the problem under study is as follows; the semi- 
implicit scheme at each time step decouples the heat (|3| and the momentum equations ^ . The time discretization 
of the heat equation is as follows: 

OiQn+l Afi n _|_ ft n ~^ 

-r— + 2u" • V0 n - u"- 1 ■ W 1 - 1 = A9 n+X , (44) 

2At y ' 

where the time derivative of the temperature field has been evaluated with a second order fixed step BDF formula. 
Once 9 n+1 , is known the velocity and pressure at time i„+i, are obtained by solving the following linear system in 
the unknown fields: 

V • u" +1 = 0, 

VP" +1 = R6 n+1 e 3 + div ( Vu " +1 + (Vu"+Y)j . (45) 
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mplicit, section 4.1 
— Semi-implicit, section 4.2 
-Scheme section 4.3 

Scheme section 4.3, variable step 



Semi-implicit, section 4.2 

Scheme section 4.3 

Scheme section 4.3, variable slep ; 




(a) Time evolution of a transitory regime 



(b) Evolutions of the error versus time 



Figure 4: Transition of an initial data to a stationary solution with T — 3.4, R — 75 for the exponential viscosity 
law Q. 



We implement this scheme by expanding the unknown fields in the equations ( 44 1 and ( 45 ) with Eq. ( 26 1 . They 



are solved at successive times t n . The discrete version of (44) is a linear system with L x M unknowns, while the 



discrete version of (45) has 3 x L x M unknowns. The decoupled nature of the procedure gives a certain speed 



advantage to this method as compared with the previous one. However, the results presented in the next section 
confirm that this method is not robust when applied to the differential algebraic equations under study. Increasing 
the order of the method or using a variable step technique does not improve the output provided by this approach. 



5 Results 

This section completes the description provided in [7] on the solutions of the convection problem in which the 
viscosity depends exponentially on temperature. By describing stationary, transitory and time dependent regimes 
for the problem under study, the consistency between the reported numerical procedures is confirmed. 

5.1 Stationary solutions and their stability 

Figure [3] displays the bifurcation diagram obtained from the analysis of a fluid layer in a finite domain with aspect 
ratio r = 3.4. Viscosity depends on the temperature according to the exponential law Q, with /i = 0.0862. The 
viscosity contrast across the fluid layer depends on R in such a way that at the instability threshold this contrast is 
around 6 • 10 2 , while at the maximum Rayleigh numbers displayed in Figure [3] it is of the order of 3.1 • 10 4 . Stable 
branches are displayed as solid lines, while unstable branches are dashed. 

The diagram displayed in Fig. [3]is obtained by branch continuation techniques, as explained in [7]. The solutions 
are obtained with the procedure described in Section 3.2. The vertical axis represents the sum \b^\ -f- \b^2 

I. These 

are coefficients obtained from the expansion of the temperature field: 

[X/21 M-l [L/2] M-l 

9(x,z,t)=J2 E b e lm (t)T m (z)cos((l-l)x)+ Z)c? m (*)T m («)siii((i-l)a;). 

1 = 1 m=0 1=2 m=0 

Most of the coefficients b lm {t) 1 c lm (t) in this expansion are approximately zero, but others are not. As a represen- 
tation of the whole spatial function, we select the sum of the significant coefficients | and \b 6 12 \- 

Table [T] confirms the convergence of the eigenvalues for a stationary solution obtained at T — 3.4, R = 78. At 
this low R number, results have two significant decimal digits for expansions L = 29 x M = 38 onwards. The results 
confirm that this is a stable solution. The presence of a zero eigenvalue is expected due to the symmetry SO (2) 
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0.1275 




(a) Time evolution of a transitory regime (b) Evolutions of the error versus time 

Figure 5: Transition of an initial data to a stationary solution with T — 3.4, R = 78 for the exponential viscosity 
law Q. 

derived from the periodic boundary conditions. Periodic boundary conditions imply that arbitrary translations of 
the solutions along the ^-coordinate must also be solutions to the system. Thus, instead of an isolated fixed point 
at the bifurcation threshold a circle of fixed points emerges [36]. The neutral direction is the direction connecting 
fixed points on this circle. Table [2] shows the convergence results obtained for R = 110. At this larger R number 
higher expansions are required; this is also expected because at large R numbers the viscosity contrast across the 
fluid layer increases. For expansions L — 55 x M — 60 onwards, a significant number of decimal digits is obtained. 

Results displayed in Figure [3] are interpreted as follows: at the bifurcation threshold the fluid undergoes a 
sub-critical bifurcation, as reported in [7J I35J. The unstable branch bifurcates at R ~ 74, below the critical 
threshold for the conductive solution, in a saddle-node bifurcation at which a stable stationary branch emerges. 
The pattern of the plume at this aspect ratio, consistent with diagram |] has wave number m = 1. The stable 
branch becomes unstable at R ~ 118 through a Hopf bifurcation. Stationary solutions are displayed at the different 
R numbers highlighted with vertical dashed lines. These results extend those reported in [7J, where the morphology 
of the plume has not been discussed. Images in Fig. [3] show the evolution of the plume with the Rayleigh number. 
As reported in [35], three idealized shapes for plumes are typically found: spout-shaped, where the tail of the 
plume is nearly as large as the head; mushroom-shaped; and balloon-shaped. Figure [3] confirms that at low R 
numbers (i? = 80) the plume is spout-shaped. At higher R = 90, 100 the plume is more rounded at the top and 
becomes closer to a balloon-shaped plume while at R = 118 the plume becomes closer to a mushroom-shaped one. 
Regarding the velocity field, Moresi and Solomatov report in [B] that from 10 4 — 10 5 viscosity contrasts a stagnant 
lid develops, and the upper part of the fluid, where the viscosity is much larger, does not move. The velocity 
fields overlapping the temperature patterns in Figure [3] confirm that at the larger viscosity contrasts obtained at 
R = 118, the velocity in the upper part of the fluid is almost null. 

The set of stationary solutions obtained with the techniques reported in Section 3, as already noted by Pla et al. 
in [7] is more comprehensive than what one would expect from mere direct time evolution simulations, because the 
latter do not prove anything about the asymptotic regime. However, as also noted by these authors, simulations of 
the evolution in time are also necessary to describe time-dependent regimes which are present at high R numbers. 
From the computational point of view, the Newton- Raphson method is more advantageous than the time evolution 
schemes, as it finds a stationary solution in less than 50 iterations, while the semi-implicit scheme needs around 
200 iterations (matrix inversions) to find the same solution beginning from the same initial data. 



5.2 Time dependent and transitory regimes 



Figure 4(a) represents the time evolution of the coefficient b%i m a transitory regime towards a fixed point. A 
rescaled Time= lOt is represented on the horizontal axis, where t is the dimensionless time. The simulation is 
produced at R = 75, which is above the instability threshold of the conductive solution, as confirmed in Figures 
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[2|and [3| Th e results displayed in Figure 4(a) are obtained with all the numerical schemes described in Section 4. 
Figure ]4(b) represents the evolution of the error versus time for the different schemes. The error is defined with 
respect to the benchmark solution obtained with the fully implicit approach. 

Despite the good performance of all schemes at low R number, at higher R numbers the semi-implicit scheme 
reported in Section 4.3 breaks. Fig 5(a) confirms this point by depicting at R = 78 the time evolution of the 
coefficient 621 m a transitory regime towards a fixed point. The schemes in Section 4.3 fail to solve the transition 
with fixed and variable time steps, which is nevertheless well determined with the implicit and semi-implicit scheme 
in Section |4.2| Fig |5(b)| reports the evolution of the error with respect the benchmark solution for the different 
schemes. 

Figure [6] confirms the existence of time-dependent convection after the Hopf bifurcation occurred at R = 118. 
Fig 6 (a) I displays the time evolution of the coefficient versus time. Snapshots of the plume in the time-dependent 



regime are shown in Figures 6(c) 6(h) In this time series, hot blobs ascend in the central part of the plume and 
these are released in the upper part of the fluid. Fig |6(b)| shows the evolution of the error for the semi-implicit 
scheme. 

As regards computational performance, the semi-implicit scheme requires higher expansions than the fully 
implicit scheme in order to achieve stability, but as it eventually requires fewer matrix inversions -both because 
it requires only one inversion at each time step and because larger step sizes are allowed- it is slightly faster than 
the fully implicit scheme. On average, for the simulations reported in this article the semi-implicit scheme requires 
80% less time than the fully implicit scheme to obtain the same results. 



6 Conclusions 

This paper addresses the numerical simulation of time-dependent solutions of a convection problem with viscosity 
strongly dependent on temperature at infinite Prandtl number. We propose a spectral method which deals with 
the primitive variables formulation. Time derivatives are evaluated by backwards differentiation formulas (BDFs), 
which are adapted to perform with variable time step. BDFs are a particular case of multistep formulas which are 
implicit. We solve the fully implicit problem and compare it with a semi- implicit method. For the problem under 
study, the proposed semi-implicit method is shown to be accurate and to have a slightly faster performance than 
the fully implicit scheme. We further show that other semi-implicit schemes, which provide a good performance in 
classical convection problems with constant viscosity and finite Pr number, do not succeed in this set-up. 

The time-dependent scheme succeeds in completing the results reported for this problem by [7j. Assisted by 
bifurcation techniques, we have gained insight into the possible stationary solutions satisfied by the basic equations. 
The morphology of the plume is described and compared with others obtained in the literature. Stable stationary 
solutions become unstable through a Hopf bifurcation, after which the time-dependent regime is solved by the 
spectral techniques proposed in this article. 

Our results confirm that the combination of direct numerical simulation, branch continuation techniques and 
linear stability analysis provide the best insight into the possible solutions of convection problems with viscosity 
strongly dependent on temperature. 
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